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Abstract 

Finite differences, as a subclass of direct methods in the calculus of variations, consist in 
discretizing the objective functional using appropriate approximations for derivatives that ap- 
pear in the problem. This article generalizes the same idea for fractional variational problems. 
We consider a minimization problem with a Lagrangian that depends on the left Riemann- 
Liouville fractional derivative. Using the Griinwald-Letnikov definition, we approximate the 
objective functional in an equispaced grid as a multi-variable function of the values of the 
unknown function on mesh points. The problem is then transformed to an ordinary static op- 
timization problem. The solution to the latter problem gives an approximation to the original 
fractional problem on mesh points. 
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^ ■ 1 Introduction 



Recently, fractional calculus, a classical branch of mathematical analysis that studies non-integer 
powers of differentiation operators, has proved a huge potential in solving complicated problems 
from science and engineering [5] [7j [H . In this framework, the fractional calculus of variations is 
a research area under strong current development. For the state of the art, we refer the reader to 
the recent book [17 , for models and numerical methods we refer to [5]. A fractional variational 
problem consists in finding the extremizer of a functional that depends on fractional derivatives 
and/or integrals subject to some boundary conditions and possibly some extra constraints. In this 
work we consider the following minimization problem: 

J[x(-)] = { L(t,x(t), a D?x(t))dt — >min, 



(1) 

x(a) = x a , x(b) = x b , 

that depends on the left Riemann-Liouville derivative, a D", which is defined as follows. 

Definition 1.1 (see, e.g., |13|). Let x(-) be an absolutely continuous function in [a,b] and < 
a < 1. The left Riemann-Liouville fractional derivative of order a, a L>", is given by 



aD?x(t) = r(l Q) - J (t-T)~ a x(r)dT, t G [a, b]. 



*This work was partially presented 16-May-2012 by Shakoor Pooseh at FDA'2012, who received a 'Best Oral 
Presentation Award'. Part of first author's Ph.D., which is carried out at the University of Aveiro under the 
Doctoral Program in Mathematics and Applications (PDMA) of Universities of Aveiro and Minho. 
Submitted 26-Aug-2012; revised 25-Jan-2013; accepted 29-Jan-2013; for publication in Computers and Mathematics 
with Applications. 
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There are several different definitions for fractional derivatives, left and right, that can be 
found in the literature and could also be included. They posses different properties: each one of 
those definitions has its own advantages and disadvantages. Under certain conditions they are, 
however, equivalent and can be used interchangeably. 

There are two major approaches in the classical theory of calculus of variations. In one hand, 
using Euler-Lagrange necessary optimality conditions, we can reduce a variational problem to the 
study of a differential equation. Hereafter, one can use either analytical or numerical methods to 
solve the differential equation and reach the solution of the original problem (see, e.g., [E]). This 
approach is referred as indirect methods in the literature. On the other hand, we can tackle the 
functional itself, directly. Direct methods are used to find the extremizer of a functional in two 
ways: Euler's finite differences and Ritz methods. In the Ritz method, we either restrict admissible 
functions to all possible linear combinations x n — y\._., ai<pi{t), with constant coefficients a; and a 
set of known base functions 4>i, or we approximate the admissible functions with such combinations. 
Using x n and its derivatives whenever needed, one can transform the functional to a multivariate 
function of unknown coefficients a^. By finite differences, however, we consider the admissible 
functions not on the class of arbitrary curves, but only on polygonal curves made upon a given 
grid on the time horizon. Using an appropriate discrete approximation of the Lagrangian, and 
substituting the integral with a sum, we can transform the main problem to the optimization of a 
function of several parameters: the values of the unknown function on mesh points (see, e.g., [5]). 

Indirect methods for fractional variational problems have a vast background in the literature 
and can be considered a well studied subject: see [TJ [2 |H [10l [T2j [T5J [20l [23] and references therein 
that study different variants of the problem and discuss a bunch of possibilities in the presence of 
fractional terms, Euler-Lagrange equations and boundary conditions. Direct methods, however, 
to the best of our knowledge, have got less interest and are not well studied. A brief introduction 
of using finite differences has been made in [22 , which can be regarded as a predecessor to what 
we call here an Euler-like direct method. A generalization of Leitmann's direct method can be 
found in [3], while |16) discusses the Ritz direct method for optimal control problems that can 
easily be reduced to a problem of the calculus of variations. 



2 Direct methods in the classical theory 

The basic idea of direct methods is to consider the variational problem as a limiting case of a certain 
extremum problem of a multi-variable function. This is then an ordinary static optimization 
problem. The solution to the latter problem can be regarded as an approximate solution to the 
original variational problem. There are three major methods: Euler's finite differences, Ritz and 
Kantorovich's methods. We are going to discuss the Euler method briefly: see, e.g., [3]. The basic 
idea of a finite differences method is that instead of considering the values of a functional 

J[ x (-)] = [ L(t,x(t),x(t))dt 

J a 

with boundary conditions x(a) = x a and x(b) = x , on arbitrary admissible curves, we only track 
the values at an n + 1 points grid, ti, i = 0, . . . , n, in the interested time interval. The functional 
J[ie(-)] is then transformed into a function ty(x(ti), xfo), ■ ■ • , x(t n -i)) of the values of the unknown 
function on mesh points. Assuming h = ti — x{ti) — Xi and Xi ~ x% one has 



J[x(-)] ~ $(x 1 ,x 2 ,...,x n -x) = (ti, 



Xq — Xa, X n — X . 



The desired values of x\, i = 1, . . . , n — 1, give the extremum to the multi- variable function ^ and 
satisfy the system 

— = 0, i = l,...,n-l. 

OXi 
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The fact that only two terms in the sum, (i — l)th and zth, depend on Xi makes it rather 
easy to find the extremum of 'J solving a system of algebraic equations. For each n, we obtain a 
polygonal line, which is an approximate solution to the original problem. It has been shown that 
passing to the limit, as h — ¥ 0, the linear system corresponding to finding the extremum of ^ is 
equivalent to the Euler-Lagrange equation for the original problem |24j . 



3 Finite differences for fractional derivatives 

In classical theory, given a derivative of certain order x^ n >, there is a finite difference approximation 
of the form 

^(*) = ^ + ^E(-i) fe 0^-^). 

k=0 

where 



*(n-l)(n-2)-.-( n-fc+l) 



kj k\ 

This is generalized to derivatives of arbitrary order and gives rise to the Grunwald-Letnikov 
fractional derivative. 

Definition 3.1 (see, e.g., [13 ). Let < a < 1. The left Grunwald-Letnikov fractional derivative 
is defined as 

1 00 / 



fe=0 



?D?x(t) = ft lim - I )*(* - kh). (3) 



Here (?) is i/ie generalization of binomial coefficients (|2|) reaZ numbers. Similarly, the right 
Grunwald-Letnikov derivative is given by 

G t L D>(f) = lim -1 £)(-!)* fjx(t + fcfc). (4) 



7^0+ /l" ^ V 

fc=0 v • 

The series in © and (0} converges absolutely and uniformly if x(-) is bounded. The definition is 
clearly affected by the non-local property of fractional derivatives. The arbitrary order derivative 
of a function at a time t depends on all values of that function in (— oo, t] and [t, oo) because of the 
infinite sum, backward and forward difference nature of the left and right derivatives, respectively. 
Since we are, usually, and specifically in this paper, dealing with closed time intervals, the following 
remark is made to make the definition clear in closed regions. 

Remark 3.2. For the above definition to be consistent, we need the values of x(t) outside the 
interval [a, b] . To overcome this difficulty, we take 



x*(t) = 



if t S [a, 6], 
if t $ [a,b). 



Thus, one can assume G ^Dfx{t) = G ^Dfx*{t) and G t L D%x(t) = G t L D%x*{t) for t 6 [a, b]. 

As we mentioned before, this definition coincides with Riemann-Liouville derivatives. The 
following proposition establishes the connection between these two definitions and that of Caputo 
derivative, another type of derivative that is believed to be more applicable in practical fields such 
as engineering and physics. 

Proposition 3.3 (see, e.g., [2I])- Let n — 1 < a < n, n e N, and x(-) € C" _1 [a, b}. Suppose also 
that x( n '(-) is integrable on [a,b]. Then the Riemann-Liouville derivative exists and coincides with 
the Grunwald-Letnikov: 



aD"x(t) = J2 Vm + - , + vUr^S / (* - r) n - l - a x^ (r)dr = G a L D?x(t). 
T(l + i-a) T(n-a)J a 



3 



Remark 3.4. For numerical purposes, we need a finite sum in §5§ and this goal is achieved by 
Remark \3.2[ Given a grid on [a, b] as a = to, t±, . . . ,t n — b, where ti — to + ih for some h > 0, we 
approximate the left Riemann-Liouville derivative as 

a D?x(ti) * E K ) x & - kh ) = : aD?x(U), (5) 
k=a 

where = (— l) fe (^) = r(-a)r{k+i) • Similarly, one can approximate the right Riemann- 

Liouville derivative by 

1 n—i 

t D?x(ti)~—J2("k)x(ti + kh). (6) 
As it is stated in this approximation is of first order, i.e., 

a D?x(U) = ^ E *(*< - kh ) + W 

k=0 

Remark 3.5. In [18], it has been shown that the implicit Euler method is unstable for certain 
fractional partial differential equations with Grilnwald-Letnikov approximations. Therefore, dis- 
cretizing fractional derivatives, shifted Griinwald-Letnikov derivatives are used. Despite the slight 
difference with respect to standard Griinwald-Letnikov derivatives, they exhibit, at least for certain 
cases, a stable performance. The shifted Griinwald-Letnikov derivative is defined by 

fc=0 

Other finite difference approximations can be found in the literature, e.g., Diethelm's backward 
finite difference formulas for fractional derivatives [3]. 



4 Euler-like direct method for fractional variational prob- 
lems 

As mentioned earlier, we consider a simple version of the fractional variational problem where 
the fractional term consists of a Riemann-Liouville derivative on a finite time interval [a, 6], The 
boundary conditions are given and we approximate the fractional derivative using the Griinwald- 
Letnikov approximation given by ([5]). In this context we discretize the functional in (JlJ using a 
simple quadrature rule on the mesh points a = to,ti, . . . ,t n — b with h — £=2, The goal is to 
find the values x\, X2, ■ ■ ■ , x n —i of the unknown function x(-) at points ti, i — 1, 2, . . . , n — 1. The 
values of xq and x n are given. Applying the quadrature rule gives 

^(•)]=E / L(U,Xi, a D?Xi)dt~J2 hL (thXha D ?Xi) 

i=l Jti-1 i = l 
and by approximating the fractional derivatives at mesh points using ([5]) we have 

n 

J[x(-)] (ti, x u a D?Xi-k) ■ (7) 

i=l 

Hereafter the procedure is the same as in classical case. The right hand side of ([7]) can be regarded 
as a function '3/ of n — 1 unknowns x = (x\, X2, ■ ■ ■ , x n -\), 

n 

*(x) = E hL (u, x h 5 t a Si_ fc ) . (8) 

i=l 
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dxi 



To find an extremum for "J, one has to solve the following system of algebraic equations: 

0, i = l,...,n-l. (9) 
Unlike the classical case, all terms in ([8]), starting from the ith term, depend on x% and we have 

(10) 



= /i— (t i ,x i , a D?x i ) + 1 ^2 — (w£) - ^— (t i+k ,x i+k , a D?x i+k ). 

k=0 a 1 X 



dxi dx 

Equating the right hand side of (fit)]) with zero one has 
dL ,, ~ n s 1 J , , dL 



dx 



(t i ,x i , a D°x l ) + j^^2i^k) q Dax (ti+k,x l+k , a D°x l+k ) = 0. 



Passing to the limit and considering the approximation formula for the right Riemann-Liouville 
derivative ([6]), we can prove the following result. 

Theorem 4.1. The Euler-like method for a fractional variational problem of the form |T]) is 
equivalent to the fractional Euler-Lagrange equation 

dL + D« dL -0 
-dx~ +tUb d^fx~- (h 

as the mesh size, h, tends to zero. 

Proof. Consider a minimizer (x\, . . . ,x n -\) of 'J, a variation function n 6 C[a, b] with 77(a) = 
77(6) = and define rji — r)(ti), for i = 0, . . . ,n. We remark that 770 = ?7„ = and that (x\ + 
£771, . . . , x n -i +e»7n-i) is a variation of (xi, . . . , x n _i), with |e| < r, for some fixed r > 0. Therefore, 
since (xi, . . . ,x n -i) is a minimizer for ^, proceeding with Taylor's expansion, we deduce that 



< ^(X! + £7/1, . . . ,X n -l + £77„_i) - . . . ,x n _i) 

"<9L r . 



; E^ 



dL 1 



k=0 



0(e), 



where here, and in what follows, 



V fe=0 / 



Since e takes any value, we can write that 



dx im+ d a D? l "'k 



'i-k 



k=Q 



= 0. 



(11) 



On the other hand, since t]q — 0, reordering the terms of the sum, it follows immediately that 



E = E * E(-*)^ + *]■ 

4=1 fc=0 i=l fe=0 ° aL>t 



Substituting this relation into equation (jlll) . we obtain 



E 7 ^ 
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Since r]i is arbitrary, for i = 1, . . . ,n — 1, we deduce that 

^H^EKl^^HO, < = 1,...,„-1. 

Let us study the case when ro goes to infinity. Let t G]a, b[ and i G {l,...,n} be such that 
tt_ i < t < tj. First observe that in such case, we also have i — > oo and n — i — >• oo. In fact, let 
i G {1, . . . , n} be such that a + (i — < t < a + i/i. Then, i < (t — a)/h + 1, which implies that 

b-t 

n — i > n- 1. 

o — a 

Therefore, 

lim = t. 

n— >-oo,z— J-oo 

Assume that there exists a function x G C[a, 6] satisfying 

Ve > OBTVVn > N : \x t -x(U)\ < e, Vi = 1, . . . , re - 1. 

As x is uniformly continuous, we have 

Ve > OBTVVn > TV : \ Xl -x{t)\ < e, Vi = 1, . . . , n - 1. 

By the continuity assumption of 3:, we deduce that 



k=Q 

For n sufficiently large (and therefore i also sufficiently large), 

dL r ., dL 

n—^oo,i—>oo 

In conclusion, 



lim — [i] = — (t,z(i), a .D£ x(t)). 



_(i,5(i), a ^(i)) + t ^_(i,^), £)^(i)) = 0. (12) 

Using the continuity condition, we prove that the fractional Eulcr-Lagrange equation (1121) holds 
for all values on the closed interval a < t < b. □ 



5 Examples 

In this section we try to solve test problems using what has been presented in Section 01 For 
the sake of simplicity, we restrict ourselves to the interval [0, 1]. To measure the accuracy of our 
method we use the maximum norm. Assume that the exact value of the function #(•), at the point 
ti, is x(ti) and it is approximated by 5$. The error is defined as 

E = max{|x(ii) — Xi\, i = 1, . . . , n — 1}. (13) 

Example 1. Consider the following minimization problem: 

a;(0) =0, x(l) = l. 
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0.2 0.4 0.6 0.8 1 

t 



Figure 1: Analytic and approximate solutions for problem (|14[) of Example [TJ 

Problem (1141) has the obvious solution of the form x(t) = t 2 due to the positivity of the Lagrangian 
and the zero value of J [£(■)]. Using the approximation 

fc=0 

for a fixed h, and following the routine discussed in Section^ we approximate problem (I14[) by 

8=1 \ fc=0 v > / 

Since the Lagrangian in this example is quadratic, system ([9|) is linear and therefore easy to solve. 
Other problems may end with a system of nonlinear equations. Simple calculations lead to 

Ax = b, (15) 

in which 





V" -1 A 2 


T.tx A i A i-i • 


Z-/£=n- 


2 AiAi— n +2 




2i=0 A t A i+l 


EL"i 2 A i 






A = 


S"=0 A t A i+2 




2-^ii—n- 






Si=0 A i A i+n-2 









where Ai = (— l) 4 /i 1-5 ( : 5 ), x = (xi, ara, • • • , a;„_i) T and b = 62, • • • , fe„_i) T wrai/i 

T7ie linear system (| X5[) is easily solved for different values of n. As illustrated in Figure [3 fey 
increasing the value of n we get better approximations. 

Although we constructed our theory for problems in the form ([T]). other operators can be 
included in the Lagrangian. Let us now move to another example that depends also on the 
first derivative and the solution is obtained via the fractional Euler— Lagrange equation. Suppose 
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that the objective functional in the fractional variational problem depends on the left Riemann- 
Liouville fractional derivative, a D" , and on the first derivative, x. If x(-) is an extremizer of such 
a problem, then it satisfies the following fractional Euler-Lagrange equation (for a proof see [T9"]): 



dL 



^7 ^ =0 



dx * b d a D?x dt V dx 
Example 2. Consider the following minimization problem: 

fi 



(16) 



J[x(-)} = / ( D t r °x(t) - x 2 {t)) dt 



ie(0)=0, sg(1) = 1. 



(17) 



In this case the Euler-Lagrange equation (|16|) gives tD® 5 ! + 2£(t) — 0. Since tD® 5 ! — — — — 



the fractional Euler-Lagrange equation turns to be an ordinary differential equation: 

1 



r(o.5) 



i(t) = 



-(1-t) 



-0.5 



2r(0.5) 

which subject to the given boundary conditions has solution 



*(*) = 



1 



2r(2.5) 



(i-ty M + l- 



i 



t + 



2r(2.5) / 2r(2.5)' 



Discretizing problem (|17p . mi/i t/ie same assumptions of Example^ ends in the linear system 



2 


-1 





•• 


• 


" 




Xi 




h 


-1 


2 


-1 


•• 


• 







X2 




b 2 





-1 


2 


-1 • • 


• 







^3 




b 3 











• • 


• -1 


2 




1 




b n -i 



(18) 



where 



0.5 

k 



i = l,2,...,n-2, 



^-5E((-D'*"C5f) 



T/ie linear system (|18[) can oe solved efficiently for any n to reach the desired accuracy. The 
analytic solution together with some approximations are shown in Figured 

Both Examples [T] and [2] end with linear systems, and their solvability is simply dependent to 
the matrix of coefficients. Now we try our method on a more complicated problem, yet analytically 
solvable with an oscillating solution. 

Example 3. Consider the minimization problem QJ with the Lagrangian 



L 



r(5.5) 



r(3.5) 



r(i.5) 



(19) 



and subject to the boundary conditions x(0) = and x(l) = 1. The functional Ldt, with 
nonnegative L, attains its minimum value for 



x{t) = 16t 5 - 20< 3 + 5t. 



(20) 
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Figure 2: Analytic and approximate solutions for problem (|17[) of Example [2j 



The appearance of the fourth power in the Lagrangian (|19p . results in a nonlinear system when we 
apply the Euler-like direct method to this problem. For j = 1, 2, . . . , n — 1 we have 



(ftSarEK 1 

i=j \ fe=0 



■)xi- k -4>(U)\ =0, (21) 



w/iere 

m r(5.5) + r(3.5) r(i.5) ' 

System (|2ip is solved for different values of n and the results are given in Figure where we 
compare the obtained approximations with the exact solution (1201) . 




Figure 3: Analytic and approximate solutions for problem of Example [3] 



6 Conclusion 

Roughly speaking, an Euler-like direct method reduces a variational problem to the solution of a 
system of algebraic equations. When the system is linear, we can freely increase the number of mesh 
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points, n, and obtain better solutions as long as the resulted matrix of coefficients is invertible. 
The method is very fast in this case and the execution time is of order 1CP 4 for Examples Q] and 
[2 It is worth, however, to keep in mind that the Griinwald-Letnikov approximation is of first 
order, 0(h), and even a large n can not result in a high precision. Actually, by increasing n, the 
solution slowly converges and in Example[2] a grid of 30 points has the same order of error, 10~ 3 , 
as a 5 points grid. The situation is completely different when the problem ends with a nonlinear 
system. In Example [H a small number of mesh points, n = 5, results in a poor solution with the 
error E = 1.4787. The Matlab built in function fsolve takes 0.0126 seconds to solve the problem. 
As one increases the number of mesh points, the solution gets closer to the analytic solution and 
the required time increases drastically. Finally, by n — 90 we have E = 0.0618 and the time is 
T = 26.355 seconds. Table Q] summarizes the results. In practice, we have no idea about the 





n 


T 


E 


Example 1 


5 


1.9668 x 10 


-4 


0.0264 




10 


2.8297 x 10 


-4 


0.0158 




30 


9.8318 x 10 


-4 


0.0065 


Example 2 


5 


2.4053 x 10 


-4 


0.0070 




10 


3.0209 x 10 


-4 


0.0035 




30 


7.3457 x 10 


-4 


0.0012 


Example 3 


5 


0.0126 




1.4787 




20 


0.2012 




0.3006 




90 


26.355 




0.0618 



Table 1: Number of mesh points, n, with corresponding run time in seconds, T, and error, E (|13[) . 

solution in advance and the worst case should be taken into account. Comparing the results of the 
three examples considered, reveals that for a typical fractional variational problem, the Euler-like 
direct method needs a large number of mesh points and most likely a long running time. 

The Euler-like direct method for fractional variational problems here presented can be improved 
in some stages. One can try different approximations for the fractional derivative that exhibit 
higher order precisions. Better quadrature rules can be applied to discretize the functional and, 
finally, we can apply more sophisticated algorithms for solving the resulting system of algebraic 
equations. Further works are needed to cover different types of fractional variational problems. 
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